Reconstructing Equilibrium Entropy and Enthalpy Profiles from Non-equilibrium 
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The Jarzynski identity can be applied to instances when a microscopic system is puUed repeat- 
edly but quickly along some coordinate, allowing the calculation of an equilibrium free energy profile 
along the pulling coordinate from a set of independent non-equilibrium trajectories. Using the for- 
malism of Wiener stochastic path integrals in which we assign temperature-dependent weights to 
Langevin trajectories, we derive exact formulae for the temperature derivatives of the free energy 
profile. This leads naturally to analytical expressions for decomposing a free energy profile into 
equilibrium entropy and internal energy profiles from non-equilibrium pulling. This decomposition 
can be done from trajectories evolved at a unique temperature without repeating the measurement 
as done in finite-difference decompositions. Three distinct analytical expressions for the entropy- 
energy decomposition are derived: using a time-dependent generalization of the weighted histogram 
analysis method, a quasi harmonic spring limit, and a Feynman-Kac formula. The three novel 
formulae of reconstructing the pair of entropy-energy profiles are exemplified by Langevin simula- 
tions of a two-dimensional model system prototypical for force-induced biomolecular conformational 
changes. Connections to single- molecule experimental means to probe the functionals needed in the 
decomposition are suggested. 

PACS numbers: Valid PACS appear here 



I. INTRODUCTION 



There exists for biomolecular systems (such as pro- 
teins or nucleic acids) substantial interest in calculating 
the equilibrium free energy difference AF between two 
states at equilibrium from repeated measurements of the 
non-equilibrium work W done on the system along irre- 
versible paths connecting the two states. One particu- 
lar avenue, part of a broader class of approaches based 
on fluctuation theorems for systems out of equilibrium 
uses Jarzynski's non-equilibrium work theorem, 
exp(— /3AF) = (exp(— /3VF)), where the angular brackets 
denote averaging over the irreversible paths that start in 
one of the equilibrium states and /3 is the inverse temper- 
ature . Hummer and Szabo have adapted the Jarzynski 
identity to show how one can calculate not just the free 
energy differences AF between two states, but also entire 
free energy profiles F{x) along some progress coordinate 
X that is being pulled under non-equilibrium. They have 
also showed that this can be done not only on the basis of 
non-equilibrium molecular dynamics (MD) simulations, 
but also using data from actual single molecule pulling 
experiments [6| . Examples involve the application of me- 
chanical forces by laser or magnetic tweezers or by atomic 
force microscopy (AFM) to drive rare transitions in single 
molecules, such as unfolding of proteins and nucleic acid 
molecules, or the dissociation of a ligand from its protein 
target. This is important because it allows the ability 
to compute free energy profiles along the direction of the 
pull in the absence of the pulling force. Such profiles can. 
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in turn, provide valuable information about the kinetics 
and mechanisms of biomolecular folding or other impor- 
tant large-scale conformational transitions. Further in- 
sight into the mechanisms of various conformational tran- 
sitions can be provided if, in addition to the free energy 
profile, one could also compute its thermodynamics com- 
ponents: entropy and enthalpy, i.e., more precisely in the 
Helmholtz representation, the entropy profile S{x) and 
the internal energy profile U{x) in F{x) — U{x) —TS{x). 
Temperature-pertubation analytical formulae to decom- 
pose free energy profiles into entropy and enthalpy exist 
for equilibrium umbrella sampling |7Hl0l| , but not for the 
non-equilibrium Hummer-Szabo approach. In previous 
work [ll| , we derived formulae for entropy-energy decom- 
positions for changes between two states in the context 
of Jarzynski's identity. The purpose of this paper is to 
extend that work and present equations in the framework 
of the Hummer and Szabo approach for the calculation 
of entire entropy and energy profiles based on an ana- 
lytical approach derived from thermodynamic identities 
involving the temperature derivative of the free energy. 
These derivatives are taken by introducing temperature 
dependance into the formula derived by Hummer and 
Szabo through the use of statistical re-weighting factors 
using stochastic path integrals techniques for Langevin 
dynamics. 



II. THEORY 

The Jarzynski identity allows for the calculation of free 
energy profiles, as described by Hummer and Szabo 0], 
using the Feynman-Kac theorem, according to the for- 
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where x denotes the puUing molecular coordinate, W — 
VF[x(t)] is the external work functional and Vp = |(x — 
vtY is the pulling potential from the AFM cantilever 
moving with velocity v. This formula can be used both 
in simulations and experiments. In the simulation, multi- 
ple trajectories are run and the work values W are com- 
puted; in the experiment, the work is integrated along 
the measured force-extension curves, which serve as one- 
dimensional trajectories. 

The fundamental origin of the ability to decompose the 
free energy profile into entropy and energy profiles stems 
from the temperature dependence of the free energy pro- 
file F{x): energy and entropy profiles can be calculated 
according to the basic thermodynamic relationships 



U{x) 
TS{x) 



d{l3F{x)) 
(dF{x) 



\ dp 



(2) 
(3) 



Because U{x) and S{x) are /3-derivatives of F{x), if one 
can reconstruct F{x) by recording trajectories at two 
nearby temperatures, one may estimate, in principle, the 
derivatives in Eqs. ([2]) and (jSj by a finite difference ap- 
proximation (if one can detect the small signal needed 
for the finite differences from noise of the data) . How- 
ever, it turns out that the derivatives can be estimated 
in closed form from trajectories recorded at a single tem- 
perature f3. This is the gist of the present contribution: 
we show how, given a sufficiently large sampling of tra- 
jectories, and assuming Langevin dynamics at /?, it is 
possible to calculate the free energy profile at any other 
value /3' in principle without approximation. This can 
be done using Eg. ([1]), in which we multiply, a la um- 
brella sampling [IJI, each functional value accumulated 
in the average estimator (■ ■ ■) by a reweighting statisti- 
cal factor, a functional $/3'[x(t)] of the trajectory x(i). 
This statistical factor will correspond to the probabil- 
ity of sampling a trajectory x{t) at the desired value 
of /3', divided by the probability of sampling the same 
trajectory at (3. In the case of Langevin dynamics, the 
trajectory probabilities can be expressed in terms of ex- 
ponentials exp(— /3yloM[x(i)]) of a stochastic action func- 
tional, ^oM[x(t)], the Onsager-Machlup action fisl IT^. 
An additional multiplicative factor of exp(— /3y(xo)) ac- 
counts for the sampling measure of the initial equilibrium 
conditions. Although the Onsager-Machlup action can 
be derived for a general form of the Langevin equation 
(even memory-dependent), let us assume, for simplicity 
an overdamped version. This is already a good approx- 
imation for molecular dynamics beyond the picosecond 
scale of biomolecular motion [l^, which is well faster 
than the dynamics of relevance for important conforma- 



tional changes. In the case of overdamped Langevin dy- 
namics (with noise S, and friction 7 obeying fluctuation- 
dissipation) , 
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(4) 



in a time-dependent force field F(x. t) = —W{x)—k{x — 
vt) including both the systematic (molecular) force de- 
rived from the non-perturbed potential V{yi) and the 
pulling harmonic force describing a virtual AFM mov- 
ing cantilever with velocity v along the pulling coordi- 
nate X, the Onsager Machlup action is the line integral 
AoM — (4m7)~-^ (m7x — F(x, t))^(it. The desired ratio 
of relative probabilities (i.e., the reweighting factor) then 
becomes 



$(/?') = exp((/3-/?')^[x(t)]). 
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where, in discrete form, 
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thermodynamic average must then be divided by the av- 
erage of the statistical factors. This produces the formula 
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where the subscript /3 on the angled brackets indicates 
that the averages are taken over a simulation or experi- 
ment run at /?. 

The /?' derivatives in Eqs. © and ([3]) can then be 

taken directly. Letting / = ^"^txfc-jf^r.Tjo)''^^^ 
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the derivatives evaluate to 
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By taking the derivatives of a, 6 and c in Eq. ()11|) with 
respect to /?' explicitly, the expression can then be eval- 
uated sA (3' — (3 providing U [x, fi) and TS{x) in terms of 
non-equilibrium trajectory averages taken at (3: 
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where ((•■•)) is evaluated as the work- weighted average 
such that ((•••)exp(-/3W^)}/(exp(-/3T^)), and the en- 
tropy profile is expressed as the difference between the 
analytical expressions for internal energy and free energy 
profiles above. Although somewhat more tediously, the 
same expressions can be obtained by explicitly showing 
the /3-dependence of all the averages (•••) in Eq. ([T|) when 
the averages are written as path integrals over all possi- 
ble trajectories with their (/3-dependent) measure writ- 
ten explicitly, followed by taking all partial /3-derivatives 
needed for Eqs. (E])-®. 

Instead of using the weighted-histogram analysis 
method in Eq. ([1]), Hummer and Szabo recently devel- 
oped a quasi-harmonic approximation formula for the 
free energy profile as a function of the molecular pulling 
coordinate 



X = vt 
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where Fp denotes the pulling force —k{x — vt). Following 
the same procedure to derive Eq. (IT2]) , we obtain the cor- 
responding quasi-harmonic approximation for the energy 
Uqiiix) and entropy 5QH(a^) profiles, as below: 
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Expressions for the energy-entropy profiles along x can 
also be derived using the thermodynamic relations and 
the Feynman-Kac theorem in the same manner that F(x) 
in our Equation ([T]) was obtained in Ref. 
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where Hf — V^(x() + mXj/2. This formula does not re- 
quire the calculations of trajectory weights. Adapting 
the weighted histogram method, the energy profile up to 
an additive constant is computed as 
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III. NUMERICAL DEMONSTRATION 

A model potential to test the method is selected for 
which free energy, energy and entropy profiles can be 
calculated analytically: 



V{x,y) =p{x) +q{x)y', 



(20) 



with p{x) = x''{x — 2)^ and q{x) = (x^ + 1), which is 
depicted in Fig. [TJa). The potential corresponds to a 




FIG. 1: (a) Two-dimensional model potential prototypical 
of biomolecular conformational changes, with conformational 
transition occurring along x across a barrier, and with a 
harmonic perpendicular degree of freedom y with decreasing 
width (decreasing entropy) as x increases. This could model 
an "unfolded" state at a; = and a "folded" state at a; = 2. 
(b) Reconstructed (-I-) and analytical (solid line) free energy 
profiles. 

bistable potential in the x direction, with minima at a; = 

and a; = 2, joined by a potential energy barrier of height 

1 at X = 1. In the y direction the potential is a harmonic 
oscillator with a force constant that varies with x, akin 
to a trough aligned along x of increasing curvature and 
with a bump in the middle. The minima at {x — 2,y = 
0) is entropically destabilized relative to the minima at 
{x = 0,y = 0) because the latter well is wider along y. 
This 2-dimensional energy is perhaps the simplest model 
for a biomolecular conformational transition (such as a 
two-state protein folder) between two states separated 
by a barrier, and with the free energy in the "important" 
direction x, (e.g., the folding coordinate) modulated by 
entropy. 

The partition function for this model as a function of 
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the X coordinate is 



Q{x) = 2m(V/3)'/'- 



-/3p(a 



(21) 



From the partition function free energy (F(x)), inter- 
nal energy [U{x)), and entropy {S{x)) profiles can be 
calculated according to the relations F — — /3^^1ogQ, 
S — — and U = F + TS. The corresponding ana- 

lytic results for these profiles are 



and 



F{x) = pix) - -log . 



U{x)^pix) + —. 



(22) 



(23) 



(24) 



To carry out a reconstruction of the energy profiles along 
the X coordinate, a harmonic pulling potential of the form 
Vp{x) = |(x — vt)"^ is superimposed on the underlying 
potential, where k is the spring constant, v is the pulling 
velocity and t is time. 

Overdamped Langevin trajectories were simulated on 
the potential described above, with parameters 7 and m 
chosen to be unity, and /? set to be 2; a number of 10^ 
trajectories were integrated for 10"* time steps of size .001 
in arbitrary units. The spring constant, k, used for the 
pulling potential was 5, and a pulling velocity of .2 was 
used so that the minima of the pulling potential would 
vary from to 2 over the course of each trajectory. When 
evaluating the quasi-harmonic approximation, k was set 
to 10, which yielded smaller deviations near the barrier. 

Free energy, energy and entropy profiles are all re- 
covered with reasonable agreement to the analytic re- 
sult. The free energy profile is exact (Fig. (Hb)). En- 
ergy profiles calculated using three different expressions 
in Eq. and ^ are shown in Fig.lUJa). Uhs{x) 

and UQ-ii{x) show some deviation at poorly sampled val- 
ues of X, around the barrier and outside of the interval 
to 2, while Ufk{x) agrees with the theoretical values 
derived analytically. Entropy profiles shown in Fig. [^fb) 
show similar behavior. C/qh and ^qh in particular have 
failed to converge for values of x over the barrier. All 
profiles have been arbitrarily aligned at a; = 0. 

The accuracy of the results shown in Fig. [1] and [H^a) 
is measured by calculating bias given as S{x) = u{x) — 

u{x), where u{x) = 1/^X)ti=i ^^^'^ Un{x) is the 

value of ri-th sample for F^isix), U-nsix), UQ}i{x), and 
Ufk{x), computed from 1000 trajectories, while u{x) 
is the exact values of the profiles derived analytically. 
N is the number of independent samples and set to 
be 1000. On the other hand, the standard deviation 
a{x) = {'^ / N Y^^^ilunix) — u{x)YY/'^ reveals the rela- 
tive precision of the values calculated with each methods. 
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FIG. 2: (a) Energy profiles reconstructed from Eq. (fT2)) . ((TSj, 
and (|19p (symbols) and analytical result (solid line) (b) Re- 
constructed (symbols) and analytical (solid line) entropy pro- 
files. 
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FIG. 3: (a) Biases 5{x) and (b) standard deviation (j{x) for 
each methods used to calculate F{x) and U{x). 



5{x) and <j{x) corresponding to four profiles are displayed 
in Fig. H 

In cases of Fy{_s{x) and L^fk, biases are negligible com- 
pared with those of Uy{_s{x) and C/qh (a;), which were com- 
puted using the trajectory weights. Likewise, the stan- 
dard deviations for the former are smaller by two orders 
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of magnitude than those of the latter in general. The 
large values of a{x) for the latter arise from the large 
variance of trajectory weights, and the increasing behav- 
ior is because of less sampling of x and the accumulated 
error in the calculations of trajectory weights. 



IV. CONCLUSIONS 

We have presented formulae, analogous to the free en- 
ergy profile formula of Hummer and Szabo, to calculate 
energy and entropy profiles. The formulae can be applied 
rigorously to driven MD simulations of biomolecules to 
providing insight into driving forces for conformational 
changes. Although not more accurate than equilibrium 
decompositions of free energy on the basis of thermody- 
namic integration or thermodynamic perturbation 
an important advantage of the decomposition of free en- 
ergy from non-equilibrium simulations over equilibrium 
simulations is that they are trivially parallel, hence many 
quick simulations can be distributed on many indepen- 
dent CPUs and run at much shorter wall-clock time. 
However, while application to simulations is straightfor- 
ward, application to single molecule pulling experiments 
is complicated by two restrictive aspects. 

The first aspect is general to reconstruction of both the 
energy-entropy profile decomposition and to the original 
free energy profile reconstruction and has to do with the 
fact that exact formulae are expressed in terms of in- 
finitely many trajectory realizations. This is in contrast 
with experimental situations, for which only finite num- 
bers of trajectory measurements are available. This con- 
vergence issue is particularly acute when the pulling coor- 
dinate is varied much more rapidly than the equilibration 
time of the coordinate. In these far out-of-equilibrium 
conditions, the averages such as those derived for inter- 
nal energy, entropy or free energy can be dominated by 
realizations that are extremely rare, and more and more 
trajectories are likely to be needed. 

The second complicating aspect has to do with the 
ability to monitor "perpendicular" degrees of freedom, 
i.e., other important degrees of freedom than the end-to- 
end extension that is being pulled. Measurement of the 
action A in Eq. ([S]) needed to reweight requires and es- 
timate of the multidimensional force F and displacement 
Ax. To work in resolving entropy profiles, more than one 
degree of freedom would need to be probed (because, for 
a one-dimensional Langevin propagation on the recon- 
structed free energy, the 1-d free energy profile is also 
the potential-of-mean-force profile (whose gradient is the 
effective force), hence the entropy profile is constant). 
While the pulling experiments typically can report forces 
and displacement along a single direction (and have been 
used in the context of fluctuation theorems to reconstruct 
free energies [T7l - [l9| ). it is possible in principle in sin- 
gle molecule experiments to assess the effect of different 
pulling geometry on the mechanical strength; this has 
been reported in two circumstances that we are aware of. 



Firstly, while most proteins studied to date involve 
pulling recombinant tandem arrays of homopolymers, 
which limits the extension geometry to that applied be- 
tween the N- and C-termini, other extension geometries 
are possible. Examples are studies of a protein, E21ip3, 
specifically labeled with a gold reactive tag at specific 
sites (20| : ubiquitin, which can form polymers between 
its C-terminus and the side chain of one of four lysine 
residues (2l| : and lysozyme, which has been linked by 
novel disulphide bonds [22] . 

Secondly, instead of (or in addition to) pulling from 
different ends of the same proteins as done in the ex- 
periments above, another experimental procedure of use 
would be to probe N- to C- terminus separation by 
pulling permutants of the protein, i.e., mutants that have 
the same sequence but different positions of the N and C 
termini; such a strategy is possible and was used to study 
the effect of protein structure on mitochondrial import 

Such type of force-spectroscopy investigations open the 
possibility to build approximate models for what ex- 
actly fiuctuates in the perpendicular manifold (i.e., which 
structural degrees of freedom collectively contribute to 
entropy as described by a one-dimensional perpendic- 
ular harmonic variable). In principle, this harmonic 
motion could then be probed spectroscopically during 
the single-molecule pulling experiment, and this strategy 
can lend the decomposition formulae presented herein 
to useful applications. Possible experimental scenarios 
might involve suitably placed two- or three-color FRET 
dyes [23, HI] or polarization-sensitive spectroscopy with 
a rigidly attached dye [2^ that can be used to gauge 
time-dependent perpendicular fluctuations concomitant 
to pulling. 

While the formulation here was based on the original 
Hummer and Szabo approach, it may be of interest to ex- 
plore the decomposition of entropy and internal energy 
in the framework of other approaches to free energy re- 
construction, such as the extensions to the fluctuation 
theorem to include reaction coordinates [27[ or the dif- 
ferential fluctuation theorem [28i] . 

V. ACKNOWLEDGMENTS 

lA acknowledges funds from an NSF CAREER award 
(CHE-0548047). 



6 



[1] Evans, D. J., E. G. D. Cohen, and G. P. Morriss. 1993. 
Probability of second law violations in shearing steady 
states. Phys. Rev. Lett. 71:2401-2404. 

[2] Evans, D. J., and D. J. Searles. 1994. Equilibrium 
microstates which generate second law violating steady 
states. Phys. Rev. E. 50:1645-1648. 

[3] Gallavotti, G., and E. G. D. Cohen. 1995. Dynamical 
ensembles in stationary states. J Stat Phys. 80:931-970. 

[4] Crooks, G. E. 1999. Entropy production fluctuation theo- 
rem and the nonequilibrium work relation for free energy 
difTerences. Phys. Rev. E. 60:2721-2726. 

[5] Jarzynski, C. 1997. Nonequilibrium equality for free en- 
ergy differences. Physical Review Letters. 78:2690-2693. 

[6] Hummer, G., and A. Szabo. 2001. Free energy recon- 
struction from nonequilibrium single-molecule pulling ex- 
periments. Proc. Natl. Acad. Sci. U.S.A. 98:3659-3661. 

[7] Postma, J. P. M., H. J. C. Berendsen, and J. R. Haak. 
1982. Thermodynamics of cavity formation in water - A 
molecular-dynamics study. Faraday Symp. Chem. Soc. 
:55-67. 

[8] Smith, D., and A. Haymet. 1993. Free-energy, en- 
tropy and internal energy of hydrophobic interactions- 
Computer simulations. J. Chem. Phys. 98:6445-6454. 
[9] Wallqvist, A., and B. J. Berne. 1995. Computer- 
simulation of hydrophobic hydration forces on stacked 
plates at short-range. J. Phys. Chem. 99:2893-2899. 

[10] Lu, N., D. A. Kofke, and T. B. Woolf. 2003. Staging is 
more important than perturbation method for computar- 
tion of enthalpy and entropy changes in complex systems. 
J. Phys. Chem. B. 107:5598-5611. 

[11] Nummela, ,]., F. Yassin, and I. Andricioaei. 2008. 
Entropy-enthalpy decomposition from nonequilibrium 
work trajectories. J. Chem. Phys. 128:024104. 

[12] Torrie, G. M., and J. P. Valleau. 1977. Non-physical sam- 
pling distributions in Monte-Carlo free-energy estimation 
- Umbrella sampling. Journal of Computational Physics. 
23:187 199. 

[13] Ousagcr, L., and S. Machlup. 1953. Fluctuations and 
irreversible processes. Physical Review. 91:1505-1512. 

[14] Machlup, S., and L. Onsager. 1953. Fluctuations and irre- 
versible process .2. systems with kinetic energy. Physical 
Review. 91:1512-1515. 

[15] Smith, J. C. 1991. Protein dynamics: comparison of sim- 
ulations with inelastic neutron scattering experiments. 
Quarterly Reviews of Btophysics. 24:227-291. 

[16] Hummer, G., and A. Szabo. 2010. Free energy profiles 
from single-molecule pulling experiments. Proc. Natl. 



Acad. Sci. U.S.A. 107:21441-21446. 

[17] Liphardt, J., S. Dumont, S. Smith, I. Tinoco, and C. Bus- 
tamante. 2002. Equilibrium information from noneequi- 
librium measurements in an experimental test of Jarzyn- 
ski's equality. Science. 296:1832 1835. 

[18] Collin, D., F. Ritort, C. Jarzynski, S. B. Smith, I. Tinoco, 
and C. Bustamante. 2005. Verification of the crooks fluc- 
tuation theorem and recovery of RNA folding free ener- 
gies. Nature. 437:231-234. 

[19] Harris, N. C, Y. Song, and C.-H. Kiang. 2007. Exper- 
imental free energy surface reconstruction from single- 
molecule force spectroscopy using Jarzynski's equality. 
Phys. Rev. Lett. 99:068101. 

[20] BrockweU, D., E. Pad, R. Zinober, G. Beddard, P. Olm- 
sted, D. Smith, R. Perham, and S. Radford. 2003. Pulhng 
geometry defines the mechanical resistance of a beta- 
sheet protein. Nature Struct. Biol. 10:731-737. 

[21] Carrion- Vazquez, M., H. Li, H. Lu, P. Marszalek, 
A. Oberhauscr, and J. Fernandez. 2003. The mechan- 
ical stability of ubiquitin is linkage dependent. Nature 
Struct. Biol. 10:738-743. 

[22] Yang, G., C. Cecconi, W. Baase, L Vetter, W. Breyer, 
J. Haack, B. Matthews, F. Dahlquist, and C. Busta- 
mante. 2000. Solid-state synthesis and mechanical un- 
folding of polymers of T4 lysozyme. Proc. Natl. Acad. 
Sci. U.S.A. 97:139-144. 

[23] Wilcox, A. J., J. Choy, C. Bustamante, and A. Ma- 
touschek. 2005. Effect of protein structure on mitochon- 
drial import. Proc. Natl. Acad. Sci. U.S.A. 102:15435- 
15440. 

[24] Xic, X. S. 2002. Single-molecule approach to dispersed 
kinetics and dynamic disorder: Probing conformational 
fluctuation and enzymatic dynamics. Journal of Chemi- 
cal Physics. 117:11024-11032. 

[25] Rasnik, L, S. A. Mckinney, and T. Ha. 2005. Surfaces and 
orientations: Much to fret about? Accounts of Chemical 
Research. 38:542-548. 

[26] Fourkas, J. T. 2001. Rapid determination of the three- 
dimensional orientation of single molecules. Opt. Lett. 
26:211-213. 

[27] Paramore, S., G. S. Ayton, and G. A. Voth. 2007. Ex- 
tending the fluctuation theorem to describe reaction co- 
ordinates. J. Chem. Phys. 126:051102. 

[28] Maragakis, P., M. Spichty, and M. Karplus. 2008. A 
differential fluctuation theorem. J. Phys. Chem. B. 
112:6168-6174. 



